Mechanism-based and computational modeling of hydrogen sulfide biogenesis inhibition: interfacial inhibition

Hydrogen sulfide (H2S) is a gaseous signaling molecule that participates in various signaling functions in health and diseases. The tetrameric cystathionine γ-lyase (CSE) contributes to H2S biogenesis and several investigations provide evidence on the pharmacological modulation of CSE as a potential target for the treatment of a multitude of conditions. D-penicillamine (D-pen) has recently been reported to selectively impede CSE-catalyzed H2S production but the molecular bases for such inhibitory effect have not been investigated. In this study, we report that D-pen follows a mixed-inhibition mechanism to inhibit both cystathionine (CST) cleavage and H2S biogenesis by human CSE. To decipher the molecular mechanisms underlying such a mixed inhibition, we performed docking and molecular dynamics (MD) simulations. Interestingly, MD analysis of CST binding reveals a likely active site configuration prior to gem-diamine intermediate formation, particularly H-bond formation between the amino group of the substrate and the O3′ of PLP. Similar analyses realized with both CST and D-pen identified three potent interfacial ligand-binding sites for D-pen and offered a rational for D-pen effect. Thus, inhibitor binding not only induces the creation of an entirely new interacting network at the vicinity of the interface between enzyme subunits, but it also exerts long range effects by propagating to the active site. Overall, our study paves the way for the design of new allosteric interfacial inhibitory compounds that will specifically modulate H2S biogenesis by cystathionine γ-lyase.

. Structure of cystathionine γ-lyase. (A) Structure of CSE showing its tetrameric organization (PDB code: 2NMP). Each subunit of CSE is color-coded differently. The PLP cofactor is shown in brown. Structural elements from the adjacent B subunit (loop  and α-helix 229-244 ) contributing residues to the active site of the A subunit are color-coded in purple. (B) Close-up of CSE active site from the A subunit. Important amino acid residues are represented in stick representation. The PLP cofactor (brown) establishes an imine bond at its C4' position with Lys212 (tan color) and is stabilized through various hydrogen bonds (represented in light blue lines) with surrounding amino acid residues, including residues from the adjacent B subunit (colored in orchid). www.nature.com/scientificreports/ catalase, has also been reported to selectively inhibit CSE but the molecular bases for such inhibitory effect have not been investigated [33][34][35] .
In this study, we report a combined kinetic, molecular docking and molecular dynamics analysis of D-penicillamine capacity to modulate the activity of CSE. Our kinetic study reveals that D-Pen inhibits both the cystathionine and cysteine cleavage reactions catalyzed by CSE through a mixed inhibition mechanism. Our molecular dynamics analysis of cystathionine binding shows that the substrate adopts a specific active configuration prior to gem-diamine intermediate formation by establishing a strong interacting network with surrounding amino acid residues and by forming H-bond contact with the O3′ atom of PLP. Our studies show that this competent state is no longer visible in the presence of D-penicillamine, which binds to the interface of two adjacent subunits. Interestingly, inhibitor binding not only induces the creation of an entirely new interacting network at the Scheme 1. Main reactions catalyzed by cystathionine γ-lyase. The α,γ-elimination of cystathionine leads to cysteine formation (top) while the α,β-elimination of cysteine produces the biological gaseous mediator H 2 S (bottom). The steps leading from the internal aldimine to the formation of the external aldimine are not shown for reasons of clarity. www.nature.com/scientificreports/ vicinity of the interface between enzyme subunits, but also exerts long range effects by propagating to the active site thus explaining its particular mode of inhibition.

Results and discussion
Purification and characterization of CSE. Human recombinant CSE was purified as described previously 24 with the following modification: after Ni-NTA affinity chromatography, the protein was either treated with TEV protease to remove the N-terminal 6 × His-tag or dialysed and then loaded onto a HiLoad 16/600 Superdex 200 column. The Q-Sepharose step was thus skipped. The purified full-length and tag-less proteins were judged to be pure ≥ 90-95% by gel electrophoresis (Fig. S1). In addition, both full-length 42 and tag-less CSE purified as a homotetramer with estimated apparent molecular masses of 152 ± 3 and 148 ± 4 kDa, respectively. Last, the absorption spectrum of purified proteins exhibits a maximum at 428 nm that is characteristic of an internal aldimine-containing CSE 24 (Fig. S1).
The goal of the present study was to estimate the effect of the thiol D-penicillamine on the kinetic parameters of CSE-catalyzed reactions. However, the influence of D-pen on the kinetic parameters associated with the α,γelimination reaction of cystathionine by human CSE could not be determined by the standard in vitro assay that uses Ellman's reagent 43 . Instead, we monitored the formation of α-ketobutyrate, another reaction product of the α,γ-cleavage of cystathionine (Scheme 1), using a discontinuous 2,4-dinitrophenylhydrazine (DNPH) activity assay. In such conditions, purified full-lenth CSE displays V max of 2.47 ± 0.05 U/mg and K M of 0.27 ± 0.02 mM with cystathionine as a substrate (Fig. S2). These values are similar to the ones determined in the classical assay (V max = 2.5-3.1 U/mg and K M = 0.28-0.40 mM) 24,42 . Tag-less CSE exhibits a V max of 2.33 ± 0.02 U/mg and K M of 0.25 ± 0.01 mM for the the α,γ-cleavage of cystathionine in the DNPH activity assay (Fig. S2). Thus, the catalytic efficiency of human CSE is not affected by removal of the N-terminal 6 × His-tag. The tag-less protein was used to perform the following studies.

Interaction of D-penicillamine with human CSE.
A previous study reported that D-pen-dependent inhibition of CSE was reversed in a concentration-dependent manner by the addition of exogenous PLP 34 . It was thus suggested that D-pen might exert its inhibitory effect on H 2 S biogenesis by interfering with the PLP cofactor. As a result, we investigated this likely interaction by fluorescence spectroscopy. Addition of D-penicillamine to human CSE results in a decrease in fluorescence emission of the PLP cofactor at 495 nm (Fig. 2), thus suggesting that D-pen binding to CSE affects the environment of this intrisic probe. Quantification of the dissociation constant between tag-less CSE and D-pen gives a K d value of 15.0 ± 4.6 mM (n = 3 ± SD), suggesting a loose interaction between both molecules.  www.nature.com/scientificreports/ Mechanism based inhibition of D-penicillamine. Then, we studied the mode of inhibition of CSE by D-penicillamine by starting with its impact on CSE-catalyzed H 2 S biogenesis from cysteine. In the lead acetate assay, D-penicillamine exhibits an IC 50 value of 2.04 ± 0.13 mM towards H 2 S production from 1 mM Cys (Fig. 3A). This result differs from the IC 50 value of 0.27 mM previously measured using the same substrate concentration in the methylene blue assay 34 . In this preceding study, D-pen was pre-incubated for 10 min with CSE prior to Cys addition, in contrast to the experiments described herein. The difference between both studies could thus come from a slow-binding inhibition process in which the enzyme-inhibitor complex is assembled on a time scale of several minutes or from an induced fit mechanism of inhibition by D-penicillamine 37 . Accordingly, we followed the inhibition potency of D-penicillamine as a function of the pre-incubating time. Nevertheless, the inhibitory capacity of D-pen is not affected by a pre-incubation period, suggesting that D-penicillamine is not a slow-binding inhibitor. The dissimilarity between both in vitro studies could therefore result in part from the differences in assay conditions and from the purity of the enzyme used. Thus, Brancaleone et al. used the methylene blue assay which displays major drawbacks 44,45 and they used a GST-CSE enzyme purified with a single chromatographic step (GST HiTrap FF affinity column) that leads to a 50-60% pure enzyme 22 . It must be noticed that Brancaleone et al. observed a significant difference in the IC 50 values exhibited by D-pen in vitro (IC 50 of 0.27 mM against recombinant CSE) and ex cellulo (IC 50 of 0.044 mM in homogenated mouse aorta samples). Whereas an induced fit inhibition mechanism could have explained such a variance, it could also result in part from the in cellulo metabolization of D-Pen into a more potent inhibitory metabolite. For example, D-pen is a substrate for D-amino acid oxidase (DAO) 46 which could thus generate dimethyl mercaptopyruvate from D-pen, the former being a substrate analog for 3-MST. Thus, either dimethyl mercaptopyruvate or D-pen itself could inhibit 3-MST, the enzyme primary responsible for H 2 S biogenesis in rodent coronary artery homogenates 47 , which is under consideration. This could also explain the results obtained ex vivo in the presence of propargylglycine since this non-specific inhibitor of PLP-dependent enzymes will deprive 3-MST of its substrate via the inhibition of cysteine aminotransferase, the enzyme involved in the biogenesis of 3-mercaptopyruvate. Last, D-pen can lead to the formation of hydrogen peroxide (H 2 O 2 ) through the oxidase activity of certain hemoproteins 35 . H 2 O 2 being involved in the regulation of rodent coronary artery blood flow 48 , the effect produced by H 2 O 2 thus generated could be added to the points already developed to partly explain the vasodilatation variation observed by Brancaleone et al. in the presence of D-pen.
Next the K i of D-pen for CSE was assessed in the H 2 S assay in the presence of varying concentration of Cys. Surprisingly, Lineweaver-Burk analysis of the data suggests that D-pen followed a mixed inhibition mechanism (Fig. S3) to impede H 2 S biogenesis, since the plots intersect outside the x-and y-axes (Fig. S4A). As a result, Michaelis-Menten plots were individually fitted with Eq. (5) to extract the dissociation constants for competitive (K ic ) and uncompetitive (K iu ) inhibition, yielding values of 1.82 ± 0.35 mM and 8.60 ± 1.22 mM for K ic and K iu , The data (n ≥ 2 ± SD) were fitted as described in "Experimental Procedures" to extract the CSE-inhibitor competitive (K ic ) and uncompetitive (K iu ) dissociation constants. www.nature.com/scientificreports/ respectively (Fig. 3B). Noteworthy, kinetic parameters for H 2 S biogenesis from Cys determined during these analyses (K M = 3.03 ± 0.44 mM and V max = 0.49 ± 0.03 U/mg, n = 3 ± SD) are like the ones previously reported 24 . Last, CSE inhibitor could harbour different inhibition constants depending on the length of the substrate, as noted with S-3-carboxypropyl-l-cysteine with regard to both cystathionine and cysteine cleavage inhibition 32 . We thus determined CSE-D-pen dissociation constants for competitive and uncompetitive inhibition during the CSE-catalyzed α,γ-elimination reaction of cystathionine (Fig. 4). D-pen also follows a mixed inhibition mechanism to hinder cystathionine cleavage, as suggested by Lineweaver-Burk analysis of the data that shows that the plots intersect outside the x-and y-axes (Fig. S4B), with values of 2.54 ± 0.65 mM and 8.55 ± 0.89 mM for K ic and K iu , respectively (Fig. 4B). D-pen is therefore slightly less capable to competitively inhibit the cystathionine cleavage reaction than the cysteine cleavage reaction. In contrast, the dissociation constant for uncompetitive inhibition is not affected by the length of the substrate. Kinetic parameters for cystathionine cleavage determined during these analyses (K M = 0.21 ± 0.02 mM and V max = 2.27 ± 0.06 U/mg, Fig. 4B) are also comparable to the ones reported above (Fig. S2B).
Altogether, our results suggest than D-pen is a mixed inhibitor of cystathionine γ-lyase (Fig. S3). Interestingly, despite its with weak potency to impair both cystathionine and cysteine cleavage activity of CSE, the inhibitory mode adopted by D-pen is quite peculiar as it differs from the ones reported for trifluoroalanine 43 , propargylglycine 36,43 and S-3-carboxypropyl-L-cysteine 32 that are suicide inhibitors. It also contrasts with the inhibitory mode detected for aminoethoxyvinyglycine, a slow and tight binding inhibitor of CSE 43 , and the mode of action of I157172 29 and aurintricarboxylic acid 30 or (R)-2-oxo-N-(prop-2-yn-1-yl)thiazolidine-4-carboxamide 31 , respectively competitive cofactor and substrate inhibitors. A minimized molecular model for cystathionine γ-lyase. Accordingly, to further investigate the molecular mechanisms responsible for the unusual mixed inhibition mechanism displayed by D-penicillamine, we performed computational modeling studies. At first, we generated a structural model of CSE based on the crystal structure of CSE (PDB code: 2 NMP). Our minimized molecular model mainly differs from the crystal structure in the following points: (1) since the quaternary structure of the enzyme consists of a dimer of dimer 36 , we choose to work with a minimalist model constituted of a single dimer comprising the A and B subunits; (2) hydrogen atoms were added such as the pyridine moiety from the cofactor and the imine bond between the Nε-lysine (Lys212) and the C4′ atom of PLP are accurately protonated; (3) the length of the aforementioned imine bond is of 1.30 Å in our model while it exhibits an unusual length of 1.96 Å in the crystal structure (PDB code: 2 NMP); (4) last, Asn161 is no longer in H-bond contact with the O3' atom of PLP.
To investigate the stability of our model, we then performed 50 ns molecular dynamics (MD) simulations and monitored the time evolution of the root mean square deviation (RMSD) and the per-residue root mean square fluctuation (RMSF) (Fig. S5). Our dimeric model is relatively stable during the time course of the simulations, and most importantly nearly all the structural elements constituting the active site fluctuate marginally during the simulations. Noteworthy, the two distinct changes in the RMSD values seen at the beginning (t = 0-0.1 ns) and at the end (t ≥ 20-30 ns) of the simulations are respectively related to the equilibration protocol and mostly to www.nature.com/scientificreports/ the movements of the N-terminal portion of the protein backbone that is usually implicated in the dimer-dimer interface, as shown in the per-residue RMSF (Fig. S5).

Computational modeling of cystathionine binding reveals a likely active site configuration prior to gem-diamine intermediate formation.
Next, we used our minimized model to perform docking and molecular dynamics (MD) simulations on cystathionine (CST) binding to the active site of CSE (Figs. 5 and S5-S9). The representative docking pose from cystathionine bound to CSE (CDOCKER energy and CDOCKER interaction energy for docked cystathionine of − 67.8 kcal/mol and − 61.1 kcal/mol, respectively) shows that essentially four structural elements (named I to IV) contribute to the substrate access channel and the active site geometry, and also provide amino acid residues for PLP and/or substrate binding (Figs. 5A and S6). www.nature.com/scientificreports/ Thus, the cohesion and opening of the substrate access channel are mainly well-maintained thanks to the hydrogen bonding network created between the side chains of amino acid residues from structural elements I and IV (*Gln49 [an asterisk denotes an amino acid residue from the B subunit] and Glu339), from structural elements I, II and III (*Ser63, *N241, *S242 and R119), and from structural elements III and IV (Arg122, Gln123 and Asp363) (Fig. 5A). Noteworthy, this network is also completed by several hydrogen bonds established between amino acid residues from the same structural element that will not be discussed in here. Besides, the PLP cofactor is involved into π-π stacking with Tyr114, its N1-pyridine-H + forms an H-bond with Glu187, and PLP is anchored by strong hydrogen bonding between its phosphate moiety and amino acid residues from the A (Gly90, Leu91, Ser209) and B (*Arg62) subunit (Fig. S6A). Additionally, the distal part of the substrate principally initiates hydrogen bond contacts at the entrance of the substrate channel with the side chains of *Arg62, *Asn241, Tyr114 that was postulated to act as a base for deprotonating the substrate 49 , Glu339 that has been proposed to govern the specificity of CSE-catalyzed reactions 50,51 , and the backbone of Thr355. Last, Asn161 and Thr355 provide hydrogen-bond interactions with the proximal α-carboxyl group of cystathionine at the bottom of the substrate binding site, as does Arg375 according to previous docking experiments performed with yeast CSE 49 (Fig. S6B). Interestingly, the H-bond contacts implicating Asn161 and Arg375 might be preserved along the course of CSE-catalyzed reactions as the carboxylate of an aminoacrylate bound to CSE (PDB code: 6NBA) forms hydrogen-bonding interactions with both amino acid residues 32 . Also, site-directed mutagenesis experiments showed that Arg375 is important for maintaining the H 2 S-producing activity of the enzyme 51 . Last, cystathionine is further stabilized by establishing van der Waals contacts with *Tyr60 and Tyr114, C-H…S and O…S interactions with the side chain of Thr355, and unconventional C-H…O and C-H…N contacts with Leu341 or Lys212, respectively.
The representative docking pose was then subjected to 50 ns molecular dynamics (MD) simulations to investigate the stability of the CSE-CST binary complex (Figs. 5B-D, S5, S7A, S8A and S9A). While CST is stable along the entire sampling interval, its binding to CSE enhances the stability of the enzyme, as shown by the reduced RMSD (∆RMSD = 1.67 Å over the last 20 ns) and RMSF (∆RMSF = 0.60 Å) values observed for the CSE-CST binary complex in comparison to the ones detected for the minimized model ( Fig. S5A and S5E). MD simulations reveal changes at the vicinity of the binding site, particularly in the initial hydrogen bonding network observed in the docking pose that maintain the unity of the active site and substrate access channel (Fig. 5B), in the interacting network of PLP that is now better stabilized thanks to its interaction with two new amino acid residues (*Tyr60 and Thr211) (Fig. S7A) as well as in the interacting network of cystathionine (Figs. 5C, D, S8A and S9A).
Importantly, the connexion between structural elements I-II, II-III and III-IV are essentially preserved despite differences with the original interacting network observed in the docking pose and regardless of structural changes in elements I-IV implicated in the substrate access channel and the active site geometry. For example, the tip of the flexible loop 44-65 (structural element I) from the B subunit moves away from its original position in part due to the breaking of the connexion between the side chains of *Gln49 and Glu339 (Fig. 5B). The later now becomes indirectly implicated in the stabilization of the substrate at the active site through H-bonding with Thr355 (Fig. S8A). In addition, the guanidinium part of Arg119 flips about 90° to form conventional H-bonds with the distal α-carboxylate of CST while still staying connected with element II through H-bond contacts between its NH2 and Nε and the side chain of *Ser242. Also, elements II and III are connected through H-bonding between *Ser242 and Gln123, itself in H-bond contact with the backbone of Arg119 (Fig. 5B). Last, element IV undergoes a rearrangement during which (1) a new β-sheet antiparallel to the existing one is formed, (2) several residues (e.g. backbones of Ser340 and Leu341, side-chains of Thr355 and Arg375) now fall within a radius of hydrogen interaction with Glu349 localized on the newly formed anti-parallel β-sheet, (3) Arg375 establishes conventional H-bonds with the proximal α-carboxylate of the substrate and hydrogen-bond contact with Asn161 (Figs. 5B, S8A and S9A).
Furthermore, cystathionine is no more H-bonded to Tyr114 that establishes H-bonding contact with *Arg62 as well as S/π and O···S interactions with the substrate (Figs. S8A and S9A). In addition, another as S/π interaction is established between *Tyr60 and the substrate (Fig. S9A). Noteworthy, the distance between the O atom from Tyr114 and the S atom from the substrate (2.93 Å) is inferior to the sum of the van der Waals radii of O (1.52 Å) and S (1.80 Å) 52 , suggesting a stabilizing interaction between these atoms. The microenvironment created around the tyrosine hydroxyl group by the interactions described above may profoundly perturb its intrinsic pK a thus rendering Tyr114 prone to participate into proton transfer during cystathionine cleavage, for instance into the protonation of the leaving group as suggested for CGL from L. plantarum 53 . Furthermore, these structural changes cause the access channel from the substrate to the active site to become slightly more constrained, as indicated by the distances between the three amino acid residues *Ser63, Ala357 and Asp363 (perimeter of 38.0 Å [docking] versus 33.6 Å [MD]) (Table S1). Also, cystathionine slides 1.3 Å deeper into the substrate access channel, which may facilitate its interaction with active site amino acid residues and shield the active site from solvent exposure. Altogether, these structural rearrangements and the creation of new interacting networks for the PLP and substrate possibly promote the organising of the active site prior to catalysis.
Strikingly, the monitoring of the average distances between the O3′ or C4′ atoms of PLP and the proximal substrate amino group unexpectedly reveals the formation of a H-bond contact between the O3′ atom and the amino group of cystathionine during MD simulations (Fig. 5C, D and Table S1). This would allow the O3′ atom to participate into substrate deprotonation though two possible mechanisms, as previously reported 54,55 . First, the O3′ atom could directly deprotonate the proximal substrate amino group as observed in L-serine hydratase from Xanthomonas oryzae pv. oryzae 54 (Fig. S10). Thus, the keto form of the O3′ atom (species A) would play a catalytic role during the transimination reaction by taking part in the proton transfer from the substrate to the lysine amino group implicated into Schiff base formation. The initial abstraction of a proton from the substrate amino group will generate the enol form of the O3′ atom (species B), thus disturbing the original internal H-bond between the O3′ atom and the N-H + from the Schiff base (species A). This would allow the Schiff linkage to www.nature.com/scientificreports/ adopt a single-bond resonance form and to easily rotate around the C4′-C4(pyridine) axis, thus creating a more electrophilic carbon at the C4′ position and facilitating the nucleophilic attack of the amine group of the substrate (species C). Second, the O3′ atom could intervene in the indirect deprotonation of the substrate as suggested by QM/MM calculations during the degradation of l-methionine by methionine γ-lyase from Pseudomonas putida 55 . In this case, the O3′ atom contributes to the water-assisted deprotonation of the substrate, which has been found more favourable and less energetic than the direct deprotonation of l-methionine by the O3′ atom to generate the nucleophilic species prior to gem diamine formation. A water molecule located 4.3 Å away from the O3′ atom in the resting state of CSE (PDB code: 2NMP) could accomplish this task admitting that such a mechanism is followed by CSE during the α,γ-elimination of cystathionine. Alternate and/or complementary functions may be attributed to the formation of the H-bond contact between the O3′ atom and the nucleophilic substrate since our modeling studies have been performed with the substrate amino group already deprotonated. The H-bond linkage could first have a self-contradictory role. On the one hand it could stabilize the nucleophile, making it less reactive and thereby increasing the intrinsic barrier for the formation of the gem diamine intermediate. On the other hand, it could better position the nucleophilic amino group in space and therefore modulate its selectivity and favour its attack towards the C4′ atom, consequently decreasing the intrinsic barrier for the formation of the gem diamine 56,57 . Second, the hydrogen bond linkage could stabilize the adducts formed between PLP and the substrate, intermediates, or the product during catalysis since its occurrence has been detected in the crystal structures of the external aldimine formed between cystathionine and PLP of L. plantarum CGL (PDB code: 6LE4) 53 and in human CSE harbouring an aminoacrylate intermediate (PDB code: 6NBA) 32 . Last, the H-bond could also modulate the electrophilic character of the C4′ atom. Thus, while the protonation state of N1-pyridine determines the absolute acidity of the Cα and C4′ of aldimines and ketimines, the protonation state of imine and O3' modulates the relative acidities of the Cα and C4′ 58,59 . Accordingly, the snapshot captured at 50 ns during MD simulations ( Fig. 5B and D) could represent an enzyme conformation needed for promoting catalysis and particularly for lowering the activation energy of the transition state leading to gem diamine intermediate formation.

Prediction of D-penicillamine binding sites and docking of the ligand.
Next, we performed computational studies to identify putative ligand binding sites within CSE and further carry out similar studies on the ternary complex CSE-substrate-ligand. The prediction of the ligand binding sites was performed with the Galaxy site program using the crystal structure of CSE (PDB code: 2NMP, A and B chains) and with the Cavity search program from Biovia Discovery Studio (DS) 2021 using our minimized molecular model as input. Molecular docking was performed with CDOCKER software implemented in DS. Best binding site were retained based on docking score (CDOCKER interaction energy and CDOCKER energy scores). Thus, amongst the seventy-one putative ligand-binding sites, eighteen of them were selected for docking studies based on their size as well as their location on the enzyme. Finally, three potential ligand binding sites in which D-pen binds with similar docking score were selected (Fig. 6). Interestingly, the 3 binding sites are all interfacial ligand binding sites. Thus, site 42 is located above the substrate access channel in contact with the flexible loop 44-65 of the B subunit (structural element I) (Figs. 6 and S11). Site 12 is located behind the active site, at the interface of α-helix 230-243 (named here α1) and α-helix 115-127 (named here α2) from the A and B subunits, respectively (Figs. 6 and 7A). Last, site 33 is located below the substrate access channel in contact with the α-helix 231-241 from the B subunit (structural element II) and the α-helix 115-129 from the A subunit (structural element III) (Figs. 6 and 8A).
Thereafter, the binding mode of D-pen was thoroughly examined. Thus, docked D-pen to site 42 (CDOCKER energy and CDOCKER interaction energy of − 35.6 kcal/mol and − 27.4 kcal/mol, respectively) establishes H-bonds with the backbone of Met354 and the side chains of *Gln49 and Glu339 (Figs. 6 and S11). D-pen also forms C-H…O interactions with the side chain of Met354 and the backbones of *Gly53 and Met354. Docked D-pen to site 12 (CDOCKER energy and CDOCKER interaction energy of − 34.1 kcal/mol and − 26.7 kcal/mol, respectively) intervenes in salt bridge formation with *Glu127 from α2, in van der Waals contacts with Arg235 and Phe238 from α1, in O…S contact with *Glu127 and in C-H⋯S interaction with Phe238, the later interaction resembling conventional hydrogen bonds 60 . Interestingly, while D-pen bound to site 12 is located more than 20 Å away from the active site, the interfacial helix α1 it interacts with establishes hydrogen bond contacts with nearby structural elements (α3 and β1) connected to the active site (Figs. 6 and 7A). Last, docked D-pen to site 33 (CDOCKER energy and CDOCKER interaction energy of − 35.7 kcal/mol and − 27.4 kcal/mol, respectively) intervenes into H-bond formation with the backbone of Gln123 and is in an H-bond interacting radius with the side chains of Gnl123 and Glu127 (Fig. S13A). In addition, D-pen establishes C-H…O connections with the side chain of *Phe238 and Gln123, van der Waals contact with *Leu239, O···S and C-H…S interactions with Glu127.
Surprisingly, even though sites 33 and 42 have already been identified in the literature, they have been designated as putative activator binding sites 28,61 . Thus, residues surrounding site 42 have been proposed to establish electrostatic interactions with the phenoxy part of isoxsuprine, a presumed allosteric activator of CSE solely based on high-throughput virtual screening performed with the protein structure in complex with propargylglycine (PDB code: 3COG) 28 . However, isoxsuprine extends from the side of site 42 to Thr355 and established electrostatic interaction with Arg119 therefore possibly disturbing the bonding of the substrate to the active site, which raises questions regarding the supposed role of isoxsuprine as an activator of CSE. Importantly, the biochemical validation of a putative isoxsuprine-CSE interaction was not carried out and MD simulations showed that isoxsuprine was not stable during the 10 ns time course of MD simulations, its binding free energy profile strongly fluctuating and becoming positive at the end of the simulations 28 . Site 33 on the other hand has been assigned as a possible binding site for the triterpenes ursolic acid and uvaol, designed putative nonpolar allosteric activators of CSE solely based on blind docking studies. Once more, the biochemical validation of the triterpenes-CSE interaction was not performed while the docking score was low (~ 3.1 kcal/mol) 61  www.nature.com/scientificreports/ Molecular mechanisms for interfacial inhibition. Next, the best docking pose of each ternary complex CSE-substrate-ligand was subjected to 50 ns MD simulations, not only to investigate their respective stability but also to analyze the structural changes imposed by D-pen binding on the active site and cystathionine positioning. Importantly, the 10-ns all-atom MD simulations of the ternary complex formed with D-pen bound to site 42 show that the protein backbone is rapidly destabilized (Fig. S5B). Although this interfacial ligand binding site is quite interesting due to its location nearby the active site and even if we cannot entirely exclude it as a likely ligand binding site, we did not consider this complex afterwards. Both the substrate and the ligand bound to sites 12 and 33 are stable during the molecular dynamics, but the ternary complexes lose the initial stability provided by substrate binding as suggested by an increase in RMSD (∆RMSD = 0.21 Å and 2.14 Å over the last 20 ns for site 12 and site 33 vs the binary complex, respectively) and RMSF (∆RMSF = 0.30 Å and 0.81 Å for site 12 and site 33 vs the binary complex, respectively) in the ternary complexes compared to the binary complex ( Fig. S5C-E). Interestingly, MD simulations induce changes in the immediate vicinity of the interfacial D-penicillamine binding sites both at sites 12 (Figs. 7B and S9B) and 33 (Figs. 8A and S9C). In site 12, D-pen now becomes stabilized through salt bridge formation with *Glu127 and Arg235, van der Waals contact with Arg235 and Leu239, O…S and π/S interactions with Phe238, unconventional C-H…O bond and conventional H-bond formation with the backbone of Leu100 or HiS119, respectively (Figs. 7B and S9B). Interestingly, new connexions are established between α1 and β2 which carries Asp187 itself H-bonded to the pyridine ring of PLP, and between α1 and the α3-β1 part which ends with amino acid residues interacting with the phosphate group of the cofactor. These results suggest that structural changes provoked at the interfacial ligand binding site could translate to the active site through these new connexions (Fig. 7B). In site 33, D-pen establishes salt bridges with Glu127 and *Arg235, and its methyl groups participate in C-H…N contact with *Arg235, alkyl/alkyl interactions with *Leu100, *Leu239 and Val124, as well as π/alkyl interactions with *Phe238 (Figs. 8A and S9C). Importantly, many perturbations recorded in the CSE-substrate-ligand ternary complexes during the 50 ns MD simulations are not observed at both site locations in the CSE-substrate binary complex, thus revealing that binding of D-pen at either site 12 or 33 induces the creation of an entirely new interacting network at the vicinity of the interface between the A and B subunits (Figs. S12 and S13). Interestingly, interface alteration has also been observed in the context of benzothiazole derivative inhibitors of triosephosphate isomerase (TIM) from Trypanosoma crizi 62 . Thus, these interfacial inhibitors produce an allosteric effect on TIM catalytic activity by producing a constrained region at subunit interface through enhancement of intersubunit H-bonds, changing the solvent exposure of catalytic residues and altering the collective dynamics of the enzyme 62 .
MD simulations also provide evidence for long range effects of inhibitor binding propagating to the active site with the initial hydrogen bonding network connecting structural elements I-IV being profoundly perturbed (Figs. 7C and 8A), as are the cofactor (Fig. S7) and substrate (Figs. S8, S9) interacting networks. Significantly, structural elements II and III are no more directly connected while structural elements I and IV are now connected through H-bonding between *Glu59 and Glu339 (site 12) (Fig. 7C) or between *Ser56 and Ser358 (site Figure 6. Location of interfacial D-penicillamine binding sites selected for MD simulations. Shown is the position of D-pen (blue sphere) interfacial binding sites. Site 42 is located above the substrate access channel in contact with the flexible loop of the B subunit that contributes amino acid residues to the active site from the A subunit. Site 12 is located behind the active site at more than 20 Å from the PLP cofactor. Site 33 is located below the substrate access channel in contact with α-helices from the A and B subunits which participate in the active site geometry. Cystathionine (CST) and the PLP cofactor are represented in green and brown spheres, respectively. www.nature.com/scientificreports/ 33) (Fig. 8A). This results among other things in the opening of the substrate access channel in comparison to the binary complex (Table S1). At the active site, the PLP cofactor is destabilized either at its phosphate group (site 12) or at its pyridine moiety (site 33) thus explaining our fluorescence spectroscopy results which show that the attachment of D-pen to the CSE disturbs the intrinsic fluorescence of the PLP probe (Fig. 2). Thus, the phosphate group of the cofactor no longer establishes hydrogen bonds with the polypeptide chains of Gly90 and Leu91 when D-pen is bound to site 12, and the pyridine ring no longer forms H-bond with Asp187 when D-pen is bound to site 33. Consequently, the PLP can no more be activated by N-protonation in the latter situation 63 . Accordingly, www.nature.com/scientificreports/ the PLP moiety loses its net positive charge that has been estimated mandatory to initiate catalysis 64 . Besides, Tyr114 falls within a radius of hydrogen interaction with O2P and O4P atoms from the PLP cofactor (Fig. S7).
In addition, the microenvironment in which the tyrosine hydroxyl group was originally found (interactions with the S atom of the substrate) no longer exists therefore rendering Tyr114 unable to participate in a putative proton transfer reaction towards the leaving group of the substrate (Fig. S8). Last, in the presence of D-pen, the substrate does not anymore slide deeper into the substrate access channel due to direct H-bonding with amino acid residues located at the entrance of the substrate access channel such as *Ser63, *Asn241 and Glu339 (site 12) (Figs. 7D and S9B) or *Glu59 and Glu339 (site 33) (Figs. 8B and S9C). As a result, the connections initially established between cystathionine and the surrounding amino acids are completely disrupted and the distances between its nucleophilic amino group and the O3′ or C4β atoms of the PLP increases significantly ( Fig. S8 and Table S1), no longer leaving the possibility for the nucleophilic group of the substrate either to establish an H-bond with the O3' atom or to attack the C4' position of the PLP in order to form the gem diamine intermediate. Altogether, our observations indicate that sites 12 and 33 are interesting target sites for the design of selective allosteric inhibitors that will become potent pharmacological modulators of CSE functions.

Conclusion
Our findings provide a comprehensive picture of the structure and dynamics of cystathionine binding to CSE, and how they are modulated by D-penicillamine binding. Our modeling studies thus reveal that cystathionine binding promotes the organizing of the active site prior to catalysis, in particular by confining the access channel from the substrate to the active site and by establishing various interactions with surrounding amino acid residues including S/π and O···S interactions with between its S atom and Tyr114 that could therefore intervene in protonation of the leaving group. Strikingly, the proximal amino group of cystathionine also forms a H-bond contact with the O3' atom of the PLP cofactor, thus suggesting that the O3′ atom could participate in proton transfer reactions and/or increase the selectivity of the nucleophilic substrate amino group and modulate the electrophilic character of the C4′ atom. Interestingly, our biochemical analyses reveal that D-penicillamine acts as a peculiar inhibitor of both cystathionine cleavage and H 2 S biogenesis by human CSE by adopting an unusual mixed inhibition mechanism. Our MD simulations disclose that D-penicillamine binds to the interface of two adjacent subunits with inhibitor binding inducing the making of an entirely new interacting network at the vicinity of the interface between enzyme subunits. Remarkably, this interface modulation propagates to the active site and promotes an inactive state of the enzyme by perturbing PLP-enzyme, substrate-enzyme, and PLP-substrate interactions, thus explaining D-pen mode of inhibition. Overall, our studies should not only help for the design of new allosteric interfacial inhibitors against H 2 S biogenesis by cystathionine γ-lyase but also facilitate the targeted virtual high-throughput screening of libraries of putative allosteric compounds. To remove the N-terminal 6xHis-tag, the CSE protein obtained after affinity purification was dialysed against buffer A, and then incubated overnight at 4 °C with His-tagged TEV protease at a molar ratio CSE:TEV of 20:1. The mixture was then passed through a 5 mL His-Trap FF column (GE Healthcare) equilibrated with buffer A. Then, the flow-through containing CSE was subjected to dialysis and gel filtration as decribed above. The plamid pKR793-TEV coding for TEV containing a 6×His-tag was a kind gift of Ludovic Pecqueur (Collège de France, France). It was transformed into BL21-CodonPlus (DE3)-RIL competent cells (Agilent Technologies) and cells were grown overnight at 37 °C in LB medium supplemented with ampicilline (100 mg/L) and chloramphenicol (15 mg/L). Then, 1 L of LB medium complemented with the same antibiotics was ensemenced at 2% with the preculture. Cells were grown to OD 600 of 0.8 at 37 °C. TEV overexpression was induced by adding 0.5 mM IPTG, and growth was continued at 25 °C overnight. Cells were harvested as decribed above and the cell pellet was treated as previously described for hCSE. Sonication was then performed for 10 min. (10 cycles; 1 min on, 3 min off). After sonication and centrifugation, TEV was purified via the batch method using Ni-NTA Superflow resin (Qiagen), as described by the manufacturer. TEV was concentrated as above, and then stored at − 80 °C. The concentration of CSE and TEV was determined using the Bradford reagent (Bio-Rad) with bovine serum albumin as a standard.
Enzymatic assay to detect cysteine production from cystathionine. Usually, the formation of cysteine from cystathionine is measured in a continuous DTNB (5,5′-dithiobis-2-nitrobenzoic acid) activity assay 24 . However, due to the presence of D-penicillamine in the reaction mixtures, we were not able to use such assay. Instead, we followed the formation of α-ketobutyrate, another reaction product formed during the α,γcleavage of cystathionine by CSE, through a discontinuous 2,4-dinitrophenylhydrazine (DNPH) activity assay. Briefly, 500 μL of the assay mixture containing varying concentrations of cystathionine (0-3 mM) and D-penicillamine (up to 20 mM) were preincubated for 5 min in 100 mM HEPES buffer at pH 7.4 and 37 °C. The reaction was initiated by adding 2-5 μg of CSE and at the desired time points (2-5 min), 100 μL aliquots of the reaction mixture were quenched by adding 10% v/v of 10% trichloroacetic acid. Samples were centrifuged (13,000 rpm, 10 min at 4 °C) and the supernatant was mixed with 20 μL of a DNPH solution (2 mM in 2 M HCl). After 10 min incubation at 37 °C, the reaction mixture was loaded on a 96-well microplate and mixed with 280 μL of 1 M NaOH. Appropriate control experiments lacking substrate or enzyme were performed in parallel. After 10 min, the absorbance at 430 nm was recorded using a BioTek PowerWave XS microplate reader. The concentration of α-ketobutyrate in the reaction mixture was calculated using a standard curve generated with known concentrations of α-ketobutyrate. The effect of D-penicillamine on CSE was determined by estimating K M app and V max app values for CSE-catalyzed α-ketobutyrate production from cystathionine (CST) a each D-penicillamine (D-pen) concentration, and using Eq. (1).
Then, the CSE-inhibitor dissociation constants K ic for competitive inhibition and K iu for uncompetitive inhibition were calculated by plotting V max app /K M app or V max app as a function of D-penicillamine concentration, respectively, and using Eqs. (2) or (3)  www.nature.com/scientificreports/ Enzymatic assay to measure H 2 S production from cysteine. Production of H 2 S from cysteine by CSE was measured spectrophotometrically by the lead acetate assay as previously described 24 . First, the effect of D-penicillamine on the activity of CSE was determined from a reaction mixture containing 1 mM cysteine, 0.4 mM lead acetate and varying concentration of D-penicillamine (up to 100 mM) in 200 mM HEPES at pH 7.4. After 3 min of preincubation of the assay mixture at 37 °C, the reaction was initiated by adding 40 μg of CSE. Then, the formation of lead sulfide resulting from the reactivity of H 2 S with lead acetate was continuously monitored at 390 nm. The half inhibitory concentration (IC 50 ) value was obtained by plotting the relative activity of CSE (A, in per cent) as a function of the concentration of D-pen and by fitting the data with Eq. (4), where ∆A is the maximal effect of D-pen on the activity of CSE and n is the Hillslope that characterizes the slope of the curve at its midpoint.
Second, the effect of D-penicillamine on the kinetic parameters of CSE-catalyzed H 2 S generation was estimated in the same activity assay. Briefly, the reaction mixture containing varying concentrations of cysteine (1-10 mM), lead acetate (0.4 mM) and D-penicillamine (0, 1, 5 or 10 mM) was preincubated for 3 min in 100 mM HEPES buffer at pH 7.4 and 37 °C. Next, 40 μg of CSE was added to the assay mixture to initiate the reaction, which was followed at 390 nm. The molar extinction coefficient for lead sulfide previously determined was then used to obtain v H2S

24
. The CSE-D-penicillamine dissociation constants for competitive and uncompetitive inhibition were calculated by plotting v H2S as a function of cysteine concentration and using Eq. (5).
Fluorescence spectroscopy. D-penicillamine binding experiments were performed by following PLP fluorescence (λ exc = 423 nm; λ em max = 495 nm) on a fluorimeter (Shimadzu) equiped with a thermostating device. Tag-less CSE (1 μM) was incubated in 100 mM HEPES at pH 7.4 and 20 °C with successive aliquots (0.5-2 μL) of D-penicillamine stock solutions prepared in the same buffer. Fluorescence emission spectra were recorded 5 min after D-pen addition. Equilibrium dissociation constant (K d ) of D-pen to CSE was determined by plotting the relative decrease in fluorescence emission at 495 nm against the D-penicillamine concentration and by using Eq. (6), whit F 0 the minimal relative fluorescence at saturating D-pen concentration, ∆F max the maximal relative decrease in fluoresence emission provoked by D-pen binding to CSE, and n the Hill coefficient that leasures the potential cooperativity in ligand binding process.
Molecular docking and dynamics studies. The crystal structure of the PLP bound human CSE tetrameric complex was downloaded from the Protein Data Bank (PDB ID: 2NMP, https:// www. pdb. org). A dimeric structure of CSE formed by chains A and B was prepared using the Prepare Protein module in Biovia Discovery Studio (DS) 2016 with the default parameters and the CHARMM force field. Briefly, all crystallographic water molecules, except HOH 512, were removed. Bond orders were assigned, hydrogen and missing atoms were added. The protonation states on protein and PLP were adjusted at pH 7.4. The imine bond between Lys212 and C4' atom of PLP was created using the Build and Edit Protein tool of DS. Finally, the protein/PLP complex 3D structure was minimized using the Adopted Basis Newton-Raphson (ABNR) algorithm.
The 3D-coordinates of l-( +)-cystathionine and D-penicillamine were generated using the Prepare Ligand module in DS with ionization states calculated at pH 7.4. Random conformations of each ligand (8 for D-penicillamin and 26 for cystathionine) were generated for docking experiments using the BEST algorithm 38 to improve the coverage of the conformational space. Docking studies were performed using CDOCKER 39 implemented in DS 2016 and all parameters were kept at their default values. Potential binding sites in the CSE protein structure were determined using the receptor cavities tool of DS and the fpocket application 40 within the Galaxy server (https:// chemi nform atics. usega laxy. eu/). Pose clustering, based on the RMSD values, was performed using the RMSD tool of DS. Best poses within each cluster, based on the CDOCKER Interaction Energy and CDOCKER Energy scores, were further refined through the In Situ Ligand Minimization protocol of DS. Following this, the binding mode of compounds was thoroughly examined using the Analysis Ligand Poses module of DS.
To assess the structural stability of the l-(+)-cystathionine-and D-penicillamine-bound CSE, 50 ns of molecular dynamics simulations were run starting from best docking models using the CHARMM36m force field and the NAMD protocol 41 implemented in DS 2021. Protein-ligand complexes were solvated in an orthorhombic box using a TIP3P water model. Periodic boundary conditions were applied with a minimum distance of 7 Å from periodic boundary and Na+ Cl− counter ions were added to neutralize the system. Solvated complexes were